Spatial Bayesian Distributed Lag Non-linear Modelling of Temperature-Related Mortality in Birmingham: Projections under RCP Scenarios
1 Introduction
The aim of this study was to analyse the association between air temperature and all-cause mortality particularly for Birmingham and to demonstrate the application of spatial Bayesian distributed lag non-linear models (SB-DLNMs) for estimating temperature-mortality relationships at fine geographical level with limited data. Using Birmingham as a case study, this framework can be replicated by other local authorities or cities, both in UK and internationally, to guide their local public health policies.
2 Methods
2.1 Data
Daily all-cause mortality data in Birmingham were extracted from the death register table that was stored in the Birmingham and Solihull Aristotle data warehouse. The dataset covered the period between 01/01/2014 and 31/12/2024 and and included the electoral ward of residence for each deceased individual.To restrict the analysis to Birmingham residents, only those associated with a valid Birmingham ward code were included in the analysis. To avoid bias in the assessment of temperature-related mortality, individuals whose underlying cause of death was COVID-19 (ICD-10 codes U07.1, U07.2) were removed to account for the inflated deaths during the pandemic period. Therefore, a total of 94,684 deaths during the study period were used for analysis..
During the same study period, the daily maximum and minimum values of air temperature on 1 km grid across UK were obtained from the HadUK-Grid database, produced by the UK Meteorological Office. The UK-wide grid data were then spatially subset to Birmingham by cropping and masking the raster layer using the dissolved polygon to exclude all cells outside Birmingham. Daily ward-level minimum and maximum air temperature were calculated using the area-weighted average of the grid cells intersecting each ward polygon. Finally, the daily ward-level mean values of air temperature were derived from the average of the daily ward-level minimum and maximum air temperature.
To capture the future effects of excess mortality under different climate change scenarios, the bias-corrected future series of daily mean air temperature on 1 km grid across UK were obtained from NERC EDS Centre for Environmental Data Analysis, available on CEDA archive [1]. This dataset contains 4 ensemble members (01, 04, 06 and 15) representing different realisations of future climate and under each member, there are four representative concentration pathways (RCPs) namely RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5 [2]. These four scenarios correspond to the pathways that describe the future greenhouse gas concentration, where RCP 2.6 represents the best-case scenario with a lower greenhouse gas emission and lower increase in future temperature by 1.6ºC. RCP 4.5 and RCP 6.0 represent intermediate level of greenhouse gas emission with mild increase in temperature by 2.5ºC and 2.9ºC respectively. RCP 8.5 represents the worst-case scenario with a higher greenhouse gas emission and higher increase in future temperature by 4.3ºC. The same procedure was used to spatially subset the projected daily mean air temperature at Birmingham ward level.
2.2 Modelling approach
This analysis adopted the extension of the Bayesian distributed lag non-linear models (DLNMs) to the spatial Bayesian DLNMs (SB-DLNMs) as demonstrated in recent methodological study [3]. For each day \(d\) and ward \(w\) in Birmingham, let \(Y_{d,w}\) denote the observed number of deaths for all causes. The association between ambient air temperature and all-cause mortality was modelled using a Poisson likelihood with a log-link function:
\[ Y_{d,w} \sim \text{Poisson}(\lambda_{d,w}) \]
\[ \log(\lambda_{d,w})=\alpha_{s(d,w)}+\sum_{c=1}^{n_c} \beta_{w,c}\, CB_{d,w,c} \]
where \(\alpha_{s(d,w)}\) denotes the stratum-specific intercepts under the use of case-crossover design, which is a methodology framework extensively used in environmental epidemiology [4]. In order to control for long-term trends, seasonality, and ward-specific time-invariant confounding, the data were stratified by ward, year, month and day of week, such that daily temperature and all-cause mortality in a given ward were compared exclusively with observations from the same ward on the same day of the week, within the same month and year. It has to be noted that strata with zero observed deaths were excluded from the model because they do not contribute to the likelihood and provide no information for parameter estimation; \(CB_{d,w,c}\) is the cross-basis matrix that captures the non-linear relationship between temperature and mortality across a range of lags. Two separate spline basis functions were specified for the temperature and lag dimensions. With regards to the exposure-response space, cubic B-splines were selected with internal knots placing at the 10th, 75th and the 90th percentiles of the temperature range. For the lag-response space, natural cubic splines were used with three internal knots placed at equally spaced values on the log scale, including an intercept. A lag window of 0-21 days was chosen to capture the delayed effect of cold and the short-term mortality displacement from the heat effect. The specification for both basis functions followed previously established practice in the literature [5].
\[ \beta_{w,c} = \beta_c + \beta^{*}_{w,c} \]
\[ \beta^{*}_{w,c}=\frac{1}{\sqrt{\tau_c}}\left(\sqrt{1-\phi_c}\, v_{w,c}+\sqrt{\phi_c}\, u_{w,c}\right) \]
Modelling for a single local authority at fine geographical level often results in extremely high variance and unreliable estimates. To mitigate this, it is assumed that neighbouring wards share similar characteristics that contribute to their underlying mortality risk. The spatial dependency was modelled within a Bayesian hierarchical framework using the Besag–York–Mollié 2 model (BYM2), by incorporating the spatial structure directly into the ward-specific cross-basis coefficients \(\beta_{w,c}\) of the cross-basis matrix \(CB\). To define spatial neighbourhood, queen contiguity was specified such that wards were considered as adjacent if they shared either a long border or just a single corner. Specifically, each coefficient \(\beta _{w,c}\) is modelled as the sum of the global fixed effect (\(\beta_{c}\)) for the entire region, and a ward-specific random effect \(\beta^{*}_{w,c}\) which accounts for the spatial dependence between neighbouring wards. The spatial neighbourhood structure under BYM2 model can be separated into scaled structured (\(u_{w,c}\) ) and unstructured (\(v_{2,c}\)) spatial components. The proportional contribution of these components is controlled by a mixing parameter \(\phi_{c}\), and the overall variability of the spatial effect is controlled by a precision parameter \(\tau_{c}\) [6].
2.3 Prior specification
Since this analysis was conducted within a Bayesian framework, prior distributions of all model parameters were specified as follows:
\[ \alpha_{s(d,w)} \sim \mathcal{N}(0,\; 10^{4}) \] \[ \beta_c \sim \mathcal{N}(0,\; 10^6) \]
\[ \Pr\!\left(\frac{1}{\sqrt{\tau_c}} > 1\right) = 0.01 \]
\[ \Pr\!\left(\phi_c < 0.5\right) = \tfrac{2}{3},\qquad0 \le \phi_c \le 1 \]
It is assumed that the strata specific intercept (\(\alpha_{s(d,w)}\)) and the global fixed effect (\(\beta_{c}\)) follow the weakly informative Gaussian priors. The prior specification for the scaled structured and unstructured spatial components of the BYM2 model uses the Penalised Complexity (PC) priors which were introduced as a interpretable and meaningful specification [7]. For the precision parameter \(\tau_{c}\) , the PC prior was set such that \(\Pr\!\left(\frac{1}{\sqrt{\tau_c}} > 1\right) = 0.01\), implying that substantial there is only 1% prior probability of large residual variation in temperature-mortality risk across wards. For the mixing parameter \(\phi_{c}\), the PC prior was set such that \(\Pr\!\left(\phi_c < 0.5\right) = \tfrac{2}{3}\), reflecting a prior specification that there is a 67% probability of less than half of the residual variation in temperature-mortality risk across wards is spatially structured. This means that unstructured, ward-specific random noise accounts for a greater proportion of the variability. These conservative priors were adopted because BYM2 model naturally shrinks towards \(\phi_{c}=0\) unless the data provide strong evidence of genuine spatial dependency [8].
2.4 Statistical analysis
The associations between temperature and mortality across all wards were estimated with a single-stage approach with the SB-DLNM. The model was fitted using the R-INLA (Integrated Nested Laplace Approximation) package which provides computational efficiency in execution times than using Markov Chained Monte Carlo method from other R packages. Given the high dimensionality of the model induced by the cross-basis matrix and ward-specific spatial effects, inference was conducted using empirical Bayes strategy instead of the full Laplace approximation.
The output of the SB-DLNMs were obtained by drawing 1,000 independent samples from the posterior distributions of the model parameters for each ward. These posterior samples were then used to derive the ward-specific temperature-mortality associations by aggregating the coefficients of the global fixed effect for the entire region and the ward-specific cross-basis coefficients. In order to calculate ward-specific relative risk (RR) curves for each ward, grids of ward-specific temperatures were constructed using the observed minimum and maximum temperatures within each ward. These grids were then used to generate prediction cross-basis matrices for each ward. These matrices were then combined with the corresponding prediction cross-basis matrices for each ward to reconstruct the cumulative temperature–mortality association across the observed temperature range. Initially, the reference point of the temperature was selected at the 90th percentile which is around 17 to 18 °C that matches the nationwide reference under UK context.
To obtain the minimum mortality temperature (MMT) and the MMT distribution, the 1,000 RR curves calculated from the previous step were re-centred by finding their minimum RR for each ward. The search for the minimum RR on each curve was constrained within the 1st and the 99th percentiles of the ward-specific temperature grid. The constraint was justified because extreme temperatures at extreme tails of the distribution are rare, excluding them can provide more stable estimates, which is in line with previous studies [9]. Once the minimum RR values were identified, the corresponding temperatures from the temperature grids were extracted and defined as the MMT. The posterior distribution of the MMT for each ward was the obtained by identifying the MMT across all posterior samples.
Once ward-specific MMTs had been identified for each posterior RR curve, excess mortality attributable to heat and cold were then calculated using the backward attributable fraction (b-AF) framework [10]. The posterior distributions of the cumulative temperature-mortality association across lags that were centred at ward-specific MMTs were reconstructed by combing the posterior samples of the cross basis coefficients. The cumulative log-relative risks were calculated by summing lag-specific contributions over the 21 lags and transformed into b-AF for each posterior draw for each ward. These b-AF(s) were the multiplied by the daily observed deaths to derive the attributable numbers (ANs) for each posterior draw for each ward. Daily temperatures that were above the MMT(s) were counted as cold-related and above as heat-related.The cold- and heat-related excess mortality were calculated by aggregating the ANs separately for cold and heat and by year. Finally, the annual cold- and heat-related excess mortality were obtained by averaging over the study period, and standardised to per 100,000 in population using the 2022 ward level population estimates.
For projections, 1,000 posterior distribution of the coefficients of the ward-specific temperature-mortality associations from the fitted model were extracted along with the ward-specific baseline daily mortality for each ward obtained by refitting the model to the observed series and averaging the resulting expected deaths. Since there ward-level population projection were not available, this analysis did not model future demographic change or adaptation, implying that the the results should be interpreted as the potential impacts of climate-driven change under a constant population and vulnerability assumption rather than predictions. The baseline impacts were calculated over 2014-2024 (using the same calender window as the study period) and extracted from each climate scenario. Future impacts for each scenario were then estimated for the 2030s, 2040s, 2050s, 2060s and 2070s decades. The procedure to obtain the future ANs and standardised excess mortality was the same as described in the previous paragraph. Given that each climate scenario projection consists of 4 ensemble members, the results were pooled posterior draw across all ensemble members.
Uncertainty was quantified using posterior median and 95% credible interval (CrI) obtained from the posterior distributions. All data processing was performed in R version 4.4.2. All statistical analyses were performed in R version 4.4.2. All data visualisations were performed in R version 4.4.2. The code is available on ….
2.5 Results
Figure 1 shows the overall ward-specific cumulative relative risk curves for an adjacent set of neighbouring wards in Birmingham. It is observed that these RR curves show the usual J-shape with varying MMT but generally occur at a relatively warmer temperatures at around the 90th percentile temperatures. Elevated mortality risks in both cold and heat temperatures are observed across wards.
Figure 2 shows the spatial distribution of the median RR at 1st and 99th percentile of the ward specific temperatures referenced to their own MMT. Billesley exhibits the highest RR at 1st percentile temperature at 2.58 (95% CrI: 1.47-4.61). In contrast, the RR values at 99th percentile are also elevated where ward Glebe Farm & Tile Cross shows the highest RR at 1.83 (95% CrI: 0.99–3.24). The cold-and heat- extreme map shows consistent elevated RR across wards, with similar RR values between adjacent neighbours, implying spatial clustering and vulnerability to lower and higher temperatures is not uniform across Birmingham.
